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The transport of charged particles in disorganised magnetic fields is an important issue which 
concerns the propagation of cosmic rays of all energies in a variety of astrophysical environments, 
such as the interplanetary, interstellar and even extra-galactic media, as well as the efficiency of 
Fermi acceleration processes. We have performed detailed numerical experiments using Monte- 
Carlo simulations of particle propagation in stochastic magnetic fields in order to measure the 
parallel and transverse spatial diffusion coefficients and the pitch angle scattering time as a function 
of rigidity and strength of the turbulent magnetic component. We confirm the extrapolation to 
high turbulence levels of the scaling predicted by the quasi-linear approximation for the scattering 
frequency and parallel diffusion coefficient at low rigidity. We show that the widely used Bohm 
diffusion coefficient does not provide a satisfactory approximation to diffusion even in the extreme 
case where the mean field vanishes. We find that diffusion also takes place for particles with Larmor 
radii larger than the coherence length of the turbulence. We argue that transverse diffusion is 
much more effective than predicted by the quasi-linear approximation, and appears compatible with 
chaotic magnetic diffusion of the field lines. We provide numerical estimates of the Kolmogorov 
length and magnetic line diffusion coefficient as a function of the level of turbulence. Finally we 
comment on applications of our results to astrophysical turbulence and the acceleration of high 
energy cosmic rays in supernovae remnants, in super-bubbles, and in jets and hot spots of powerful 
radio-galaxies. 

PACS numbers: 98.70.Sa, 95.85.Ry, 52.25.Gj, 05.20.Dd, 05.45.Jn, 05.45.Pq 



I. INTRODUCTION 



The knowledge of the transport properties of charged particles in turbulent magnetized plasmas is a long-standing 
problem, which bears directly on many astrophysical issues, such as the penetration of low-energy cosmic rays in the 
heliosphere [Q , the propagation and escape of galactic cosmic rays in and out of the interstellar mag netic field |,|, |, 
or even the efficiency of Fermi acceleration mechanisms, in particular at shocks The diffusion coefficient transverse 
to the mean component of the magnetic field plays a particularly important role in these issues, but to date, there is 
no satisfactory description of perpendicular transport. Some studies have built upon or tried to extend the results of 
the "quasi- linear theory" jsj, whose validity is limited to very low level turbulence, i.e. a turbulent component much 
weaker than the uniform magnetic field, and which calculates the transport coefficients by statistical averages of the 
displacements perturbed to first order in the inhomogenous field. Other studies have appealed to phenomenological 
approximations such as the Bohm estimate for the diffusion coefficient D ^ tlw, which corresponds to the assumption 
that the mean free path for scattering D/v of a particle of velocity v, is given by the Larmor radius tl. This 
approximation originates from laboratory experiments which led D. Bohm to the empirical formula -Db — 0.06eT/_B for 
a plasma with temperature T. A theoretical derivation of this formula was proposed later by Taylor and McNamara 
and then extended to relativistic particles 0, but no theory of Bohm diffusion (relativistic or not) in magnetic 
irregularities has been derived stricto-sensu so far. Therefore it appears that important physical and astrophysical 
issues are yet to be answered: 

• How do the transport properties change when the level of magnetic turbulence is increased? What are the 
transport properties when the mean field vanishes? Notably, what is the relevance of the Bohm scaling? 
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• Even for low level turbulence, transverse space diffusion is not well known. It nevertheless plays a crucial role 
in the confinement of cosmic rays in galaxies or other extragalactic objects (notably radio-galaxies jets). Its 
magnitude is also of direct relevance to the performance of Fermi acceleration at perpendicular shocks. 

• Do subdiffusive and more generally anomalous diffusion regimes exist? If yes, they are also of importance for 
Fermi acceleration. 

In order to shed light on these issues, we have performed extensive numerical experiments to determine the pitch 
angle scattering rate, and the parallel and perpendicular spatial diffusion coefficients for a wide range of rigidities and 
turbulence levels. Our experiments are conducted by Monte-Carlo simulations in which we follow the propagation of 
a relativistic test particle in a stochastic magnetic field constructed from three-dimensional Kolmogorov turbulence, 
and calculate the diffusion coefficients from the statistical correlations along the trajectory. Our study is similar to 
the recent work of Giacalone and Jokipii |^ in which the spatial diffusion coefhcients in two- and three-dimensional 
magnetostatic turbulence were measured using Monte-Carlo simulations for various turbulence levels and rigidities. 
Our study is however more extensive than that of Ref. |^ . In particular we measure the diffusion coefficients in a 
broader range of rigidities, by studying the diffusion of particles with Larmor radii larger than the coherence length, 
and in a broader range of turbulence levels, by going up to pure turbulence in which there is no uniform component of 
the magnetic field. In contrast, Ref. studies the case of lower-energy particles and smaller turbulence levels, with a 
turbulent magnetic field never exceeding the uniform component in strength. We also study in detail the pitch angle 
scattering rate, which is of central interest in applications to shock acceleration processes, and study in more detail the 
issue of transverse diffusion and its relation to the chaotic wandering of field lines. Finally we will repeatedly compare 
our results to Ref. ||^ where there is overlap, which is important since these numerical experiments are delicate. 

Among our results, we confirm the extrapolation to high turbulence levels of the scalings predicted by the quasi- 
linear theory for the scattering rate and the parallel diffusion coefficient at low enough rigidity. The perpendicular 
diffusion coefficient is shown to follow a law which is quite different from the predictions of the quasi-linear theory at 
low rigidities. We argue that its behavior is compatible with chaotic wandering and diffusion of the magnetic field 
lines to which particles are "attached" . In particular, we demonstrate the chaotic behavior of the magnetic field lines 
and calculate the associated Kolmogorov length and diffusion coefficient in terms of the turbulence level. We also 
show that the Bohm diffusion coefficient only holds in a limited range of rigidities O.I ^ p 1 for pure turbulence, 
and does not exist when the mean field is non- vanishing. In this latter case the Bohm value for the coefficient is 
only obtained at maximum pitch angle scattering, i.e. for particulcs with Larmor radius of order of the coherence 
scale. Finally, we also found that diffusion operates even for particles whose Larmor radii is larger than the coherence 
length, as far as we have searched in rigidity (1.5 decade). On these scales, the scattering rate decays as expected, 
albeit moderately, as the power ~ —7/3 of the rigidity. 

Our study is conducted with the following main simplifying assumptions: 

• The magnetic field is composed of a mean homogeneous field Bq and an inhomogeneous component B: B = 
Bo-hB(x). 

• The magnetic disturbances are considered to be static. This assumption is well justified as the waves propagate 
with velocities of the order of the Alfven velocity va, smaller than the velocity of particles ^ c (we consider 
relativistic particles), and the electric force is thus smaller than the magnetic force by a factor va/c. The first 
correction to the theory is the celebrated second order Fermi process which can be described by diffusion in 
momentum space, with diffusion coefficient T{p) ~ Vs,iP'v\l (? ^ with p particle momentum, and fg, the angular 
scattering frequency, is an outcome of our study. 

• The magnetic perturbations are distributed according to isotropic turbulence, whose power spectrum is written 
in terms of Fourier momentum k as: {B^k^) cx for fcmin < k < fcmax, zero otherwise, and (B(fc)) — 0, i.e. 
random phases. The exponent /3 characterizes the properties of the turbulence, and we will concentrate on the 
case /9 = 5/3 in our numerical applications, which describes Kolmogorov turbulence. The smallest turbulence 
wavenumber is related to the maximum scale of the turbulence Lmax I'io,: fcmin — 27r/L,„ax- This largest scale 
also corresponds to the correlation length of the magnetic field to within a factor of order unity [see Eq. (p^)] . 

Our notations are as follows. The quantities we will be interested in are the scattering rate or scattering time 
Ts = l/t's, defined as the correlation time of the pitch angle, the spatial diffusion coefficient along the mean field 
Z)|| and the transverse spatial diffusion coefficient D±. These coefficients are evaluated in terms of turbulence level 
■q = (B^)/(i?^) = (B^)/[i3Q -|- (B^)], and rigidity p = 27rrL/imax = ''Lfcmin- For convenience, the Larmor radius tl is 
defined with respect to the total magnetic field: tl = e/ZeB for a particle with energy e and charge Ze. The Larmor 
pulsation of a particle of energy e is defined, for convenience, as wl = ZeBc/e, and the Larmor time = ('^l)""'^- We 
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define the scattering function as g{p, rf) = Vs/^'l = ^l/ts- When useful, we will denote by the Larmor pulsation in 
the mean field. 

The paper is organised as follows. In Section 2, we recall the relation between spatial diffusion and the scattering 
off magnetic disturbances and present the numerical method. In Section 3, we present our numerical results and 
discuss the issue of transverse diffusion and the measurement of magnetic chaos characteristics. A discussion with 
direct astrophysical consequences of our results is given in Section 5, and conclusions are offered in Section 6. Finally, 
in Appendix A, we propose a theoretical interpretation of the regimes of diffusion observed, and in particular of the 
existence of diffusion for Larmor radii larger than the maximum scale of turbulence. 



II. MOMENTUM SCATTERING AND SPATIAL DIFFUSION 



High energy particles interact with cosmic matter mostly through scattering on the magnetic field which is more or 
less frozen in the medium. The interaction is elastic in the frame of a magnetic disturbance and it can be considered 
as elastic in the plasma rest frame to lowest order in v\/c, if the disturbance propagates at small enough velocity 
va «C c. With respect to a given direction, chosen as that of the uniform component of the magnetic field if this 
latter is non-vanishing, the pitch angle of the particle changes almost randomly if the magnetic field is sufficiently 
disorganized (this will be made more precise further on). Therefore the position of the particle changes according to 
a random walk on a time scale which is longer than the coherence time of the pitch angle cosine, and pitch angle 
scattering is thus responsible for the diffusion of particles. However it is generally believed that transverse diffusion 
may also occur through wandering of the magnetic field lines. In this picture, the transverse velocity of the particle 
changes through resonant diffusion as before while the guiding center of the approximate helical motion wanders 
with the magnetic field line to which it is attached, and performs a random walk in the transverse direction. These 
notions will be quantified in the forthcoming sections. Our main objective here is indeed to quantify these various 
contributions to the process of diffusion. 



A. Definitions — Scattering time and diffusion coefficients 

We define the pitch angle a with respect to the mean field direction when it exists, otherwise the direction can be 
arbitrarily chosen. The convenient random function is the pitch angle cosine: fi{t) = cos (a), and a is a function of 
time. Since we assume a static spectrum of magnetic perturbations, the auto-correlation function of will become 
stationary in the large time limit. It can then be defined as: 



C(r)^(Mi + r)Mt)>/(MO'> (1) 

where the average can be performed in three different ways. In the original quasi-linear theory, this average is taken 
over the phases of the magnetic disturbances. In the theory of chaos, the average is performed over the phase space 
subset of chaotic motions. In practice, and this is what we will use in the numerical experiment, we assume ergodicity 
and make temporal average. Our procedure of calculating averages is explicitcd further below. 
The scattering time Tg can then be defined as the coherence time of the pitch angle cosine: 

drC(T) . (2) 







In particular, if the auto-correlation function falls off exponentially C(t) — exp {—t/T) then Tg = T. Turning to the 
spatial diffusion coefficient D\\ , let x\\ be the coordinate of a particle along the mean field direction. Then dccy = v^(t)dt 
with a constant velocity v (in our case v = c), since energy is conserved. Consider now a random variation Axy of xy 
during the time interval At supposed to be larger than the scattering time Tg. One has (Aa;||) ~ 0, and 

pt+At pt+At 

{Axh=v^ dh di2 (/i(ii)M^2)). (3) 



t Jt 



Beyond the scattering time Tg, if the stationary random process fi{t) explores uniformly the interval (—1, +1), the 
space diffusion coefficient parallel to the mean field stems straightforwardly from its definition, Eq. and Eq. (|^) : 

^11 = = T;v^rs ■ (4) 

" 2Ai 3 ' 
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Here as well the average can be made according to one of the three ways explained above. The main goal of the 
computation is then to determine the dependence of Tg, or the scattering function g, in terms of the rigidity p and 
the turbulence level rj. The theoretical result is known in the regime of weak turbulence [||, if the correlation time 
of the force suffered by the particle is much smaller than the scattering time Tg. To make it more precise, particles 
undergo resonances with the MHD modes such that fcyw/i ± = 0. The correlation time is related to the width of 
the resonance in the mode spectrum, such that : 

Afcii 

= A{k\\vpL±uji,) = v\p.\Ak\\ =uji,—^ (5) 

I 

where Afcy denotes the spectrum width, in the parallel direction. Since t'^^ ^ V^l^ ''c ^ ''"s is equivalent to ^ 
Afcjj/fc|j. In this case the memory of the initial pitch angle can even be kept and the scattering function g ^ 77(p|/i|)'^~^. 
However diffusion coefHcients calculated on timescales larger than Ts must be averaged over fi. 

Due to rotation invariance around the mean field direction, there is a single transverse diffusion coefficient (when 
diffusion occurs), given by: 



_ (A 
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where Aa;_L denotes the displacement perpendicular to the mean field during the time interval At. In weak turbulence 
theory (ji <ti 1), the gyro-phase of the particle is only weakly perturbed by the disorganized component of the field 
and ip uji,, where the gyro-pulsation wl is determined with respect to the mean field. The transverse velocity can 
be approximated by: 

v_L ~ w sin a{t) [ei cos -0 — sign{q)e2 sin t/"] , (7) 

where q denotes the charge of the particle, and we implicitly assumed the mean field to lie along the direction 63. 
The pitch angle sine sina(t) varies on the timescale Tg, which is much longer than the Larmor time in the weak 
turbulence regime. The time correlation function of the pitch angle is obviously the same as that of the cosines since 
(cos(ai -|- 02)) = 0, hence (sinai sin 02) = (cosai cos 02). Therefore the transverse diffusion coefficient reads: 



00 



Di_ = \v'' I dTC(T)cosKr) (8) 
-J Jo 

Assuming that the correlation function C(t) decays exponentially on the characteristic time Tg, one finally obtains 
a result similar to the so-called classical diffusion that reads 



This transverse diffusion based on pitch angle scattering only leads to the ratio 

Di_ 1 



^11 l + (A||/rL 
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(10) 



where Ay = SDy/f is the mean free path of a particle along the mean magnetic field. This relation can also be 
obtained by treating the magnetic disturbances as hard sphere scattering centers with weak or strong turbulence. It 
is also a result of the study of Ref . |^ , which estimate phenomenological diffusion coefficients by using well-motivated 
assumptions for the velocity auto-correlation functions of the particle orbit. Finally, since (cjlTs)^ ^ 1 in the weak 
turbulence regime, one expects D±^ <ti when ry ^ 1. However the transverse diffusion may turn out to be larger 
than predicted by quasi-linear theory, even for moderate turbulence. In particular note that Eq. (||) for the parallel 
diffusion coefficient rests on the sole assumption that C(r) vanishes on timescales longer than Tg, while the quasi-linear 
result for the transverse diffusion coefficient, Eq. (^, assumes that the particle orbit is only weakly perturbed and the 
timescale of variation of the pitch angle is much longer than the Larmor time, i.e., that the level of turbulence rj <^ 1. 
We refer to this result as a prediction of quasi-linear theory; it neglects the diffusion of the guide center carrying field 
line and the associated process of chaotic magnetic diffusion which has been analysed by Jokipii and Parker pX[ | , and 
to which we will come back in the following Section. Finally in all cases one should obtain D± D|| when 77 — s- 1, 
since the mean field vanishes in this limit and there is no prefered direction anymore. 



B. Numerical simulations 



In order to evaluate the transport coefficients, we follow the propagation of particles in stochastic magnetic fields 
by integrating the standard equation of motion (Lorentz force) , and measure the statistical quantities of interest to 
us, namely = 1/ts, and Dj_, using the estimators defined respectively in Eqs. (^,(||),(^). Strictly speaking, 
the averages contained in these expressions should be taken over the phases of the magnetic inhomogeneities. In 
practice however, one may as well take these averages as follows. For a given At [using the notations of Eqs. 
a time t is picked at random over the trajectory, and the correlation between the positions at times t and t + Ai is 
recorded; this operation is repeated many times and the average is kept. This latter is then further averaged over a 
population of particles with random initial positions, and then over an ensemble of magnetic field realizations with 
random phases. In practice, we propagate 20 — 50 particles, measure the correlations at 5000 — 10000 different times 
along the trajectory of each particle, and use a few magnetic realizations. This procedure allows to reach a sufficiently 
high signal-to-noise ratio in the simulation for a moderate computer time, as indeed setting up the magnetic field and 
propagating a particle is much more costly than taking averages along the trajectory. 

In principle one could as well take the average {Ax'^)/At as the variance of the displacement at time At over a 
population of particles originally concentrated at the origin, as in Ref . Q . However this method requires to follow 
the trajectory of a large number of particles ^ 10'^ in order to achieve a reasonable signal-to-noise ratio. The method 
we employ, which measures the correlations along the trajectory of each particle, before averaging over a population 
of particles, is less costly in computer time (but requires much more memory). Nevertheless, we also checked (and 
found) that the method which measures the variance of the displacement gave results in agreement with our method 
within the error attached to the small number of particles propagated. 

The magnetic field can be constructed in two different ways which both present pros and cons. The first method 
uses Fast-Fourier Transform (FFT) algorithms to set up the magnetic field on a discrete grid in configuration space, 
starting from the magnetic field defined through its power spectrum in Fourier space, i.e.: 



B(x) = K e(n)A{n) exp 

n 

In this equation, n is the tri-dimensional wavenumber vector, with integer coordinates taking values between 1 
and fcniax/2fcmin, e(n) is a unit vector orthogonal to n (this ensures VB = 0), ^(n) is the amplitude of the field 
component, and is defined such that {A{n)) — and {A{n)A*{n)) = k~^~'^, where the average concerns the phases of 
the magnetic field. Finally, k is a numerical prefactor which ensures the correct normalization of the inhomogenous 
magnetic component with respect to the mean field, by using the following ergodic approximation to (i?^): 



2inn.x 



(11) 



(i?^)^! I dxB2(x), (12) 

and as before r] = {B'^)/{Bq + (-B^)). In practice, the field components are calculated at each vertex x.; of a discrete 
grid in configuration space beforehand. The boundary conditions are periodic with period imax, and the fundamental 
cubic cell size is imax/-^g, where Ng represents the number of wavenumber modes along one direction. One thus has: 
fcmax/fcmin = imax/imin = whcrc the factor 2 comcs from the fact that one must consider both negative and 

positive fc-modes to respect the hermiticity of B(k). In our simulation, we typically use Ng = 256 and in some cases 
Ng = 512 which gives us a dynamic range of two orders of magnitude. 

During the propagation of particles, it is of course necessary to know the magnetic field at any point x for the 
integration of the equations of motion. Our numerical code calculates the value B(x) either by tri- linear interpolation 
between the known values of the field components on the 8 vertices of the cell to which x belongs, or by taking the 
value of B at the vertex closest to x, which amounts to assuming a constant B in cells of size imax/-^g centered on 
each vertex. While the former method does not respect VB = 0, the latter implies a discontinuous magnetic field on 
each cubic cell face. We will show in the following that the results obtained by these methods differ only when scales 
smaller than the cell size are concerned, as expected. 

A second algorithm for computing the magnetic field has been proposed by Giacalone & Jokipii ||] and calculates 
the magnetic field as a sum over plane wave modes (we will refer to this algorithm as GJ). The expression defining 
B(x) is very similar to Eq. (^l]) above, except that n needs not have integral coordinates anymore, as Fast-Fourier 
Transform methods are not used. Indeed, one does not calculate the field on a discrete grid beforehand, but its values 
are calculated where and when needed during the propagation directly from the sum over plane waves. Also the sum 
is not tri-dimensional, but one-dimensional; the wavenumbers directions are drawn at random, and the amplitude 
A{k) oc to account for phase space volume. In practice, it is convenient to have logarithmic spacing of the 
fc— modes between fci„i„ and fcmax- 
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FIG. 1: Self-correlation function of the pitch angle cosine shown as a function of time r (in units of Larmor time ti, = 
for various rigidity p = 0.072, 0.12, 0.19, 0.32, 0.52, 0.85 and for rj = 0.1. 

One main advantage of the GJ method is that there is no restriction in dynamic range due to memory usage, and 
consequently fcmax/fcmin can be as large as required. However, one is limited in terms of computer usage time since it 
is expensive to perform the sum over the wavenumber modes at each point of the trajectory if the number of modes 
iVptu becomes significant. In practice, Np^ = 500 is a strict upper limit for our applications and even with 
Npw = 200 the calculation is already much slower than a similar calculation with the above FFT algorithm. 

The number of modes is important as it controls the efficiency of diffusion, since pitch angle scattering proceeds 
mainly through resonance of the particle momentum on the magnetic field modes. In quasi-linear theory the resonance 
condition reads pfik\\ = ibfcmin, where fcy is the component of the wavenumber along the mean magnetic field direction. 
The FFT and GJ share a similar number of resonance modes in this limit 77 ^ 1. However for each resonant fcy the 
FFT algorithm has Ng ~ lO"' — 10^ transverse components to be compared with one for the GJ algorithm. One thus 
expects that at higher turbulence levels, diffusion should be more effective in the FFT algorithm due to the much 
larger total number of modes than in the GJ algorithm. Furthermore, in order to preserve a correct spacing of modes 
in the GJ algorithm, one cannot indefinitely increase the dynamic range fcmax/fcmin since Npw is fixed for practical 
reasons, i.e. computer time. 

However the FFT algorithm suffers from other limitations (apart from the limitation in dynamic range): the 
interpolation of B on scales smaller than the cell size, and the periodicity on the scale Lmax- These limitations are not 
present in the GJ algorithm, and imply that the results of the FFT method obtained for Larmor radii much smaller 
than the cell size, i.e. p <C ^/Ng, or much bigger than the periodicity scale, p 3> 1, cannot be trusted, since these 
regimes are likely to be dominated by systematic effects related to the discreteness or to the periodicity. 

Overall both methods appear complementary to each other, and we use them in turn to compare and discuss the 
robustness of our numerical results with respect to the assumptions made. 

III. RESULTS AND DISCUSSION 

A. Pitch angle scattering and parallel diffusion 

The first numerical investigation to perform is the self-correlation function of the pitch angle cosine. The behavior 
of this function is shown in Fig. (|l|) vs time interval r for various levels of turbulence 77. Two bumps are observed at 
one and two Larmor periods. These bumps are observed as long as the regular magnetic field Bo exists. Since the 
decorrelation times for rj < 1 are larger than one Larmor period, the Larmor motions are not completely disorganized 

and contribute to the correlation function with some harmonics generated by nonlinearities. The inflexion of the 

2 

function indicates that it behaves as e~"''^ as r — > and then decreases exponentially in e~'^'^ as r — s- cx). Thus the 
determination of Tg and i^g by a numerical integration of the correlation function is accurate. 

In Fig. (^ we show the scattering frequency g{ri, p) — Vs/uj-l, which is the main quantity of interest for evaluating 
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FIG. 2: The scattering function g{ri,p) = Vs/(jJl as a function of rigidity p. The symbols correspond to the measurements 
made through our Monte-Carlo experiments and correspond to various turbulence levels, as indicated. These results have 
been obtained using the FFT numerical method (see text), except for the filled circles which correspond to the GJ algorithm. 
The vertical dashed lines indicate the range of validity of our FFT algorithm {i.e. all symbols except crosses), delimited 
by /Omin = femin/fcmax, aud Pmax = Ztt, which correspoud respectively to Larmor radii t-l = imax/'rA'^g (V''" cell size) and 
tl = imax- The simulation for rj — 1 shown by crosses has been obtained with a much larger dynamic range than the others, 
i.e. fcniax/fcmin = iC. Finally, the dotted lines correspond to power law approximations with slopes 2/3 and —4/3. See text 
for comments. 



the transport coefficients. This figure shows several interesting features which deserve further comments. First of all, 
one finds that both methods for calculating the magnetic field, i.e. FFT and GJ, agree well within the range of validity 
of the former method, namely for pmin ^ p ^ Pmax, where pmin = fcmin/fcmax, and pmax = 27r. These two limiting 
rigidities correspond to Larmor radii of order of the cell size and of the maximum scale of turbulence respectively, and 



result from the discreteness and periodicity of the magnetic field grid, as explained in Section [IB. One finds that the 
scattering function behaves as a power-law with different slopes depending on the rigidity and turbulence level. For 
P ^ Pmin, it must be emphasized that the results cannot be trusted for the FFT results, i.e. all symbols except filled 
circles, and the change of slope may be artificial. For 77 < 1 and p < 1, it appears that g{r], p) oc rjp^^^, in accordance 
with the quasi-linear prediction since 2/3 = /3 — 1. 

For p > 1, one finds g{ri,p) oc r/p~*/^, an unexpected result, since the resonance conditions cannot be satistified at 
these high rigidities, and the quasi- linear theory thus predicts a sudden drop of the scattering frequency at p > 1. 
In Appendix A, we provide a first theoretical explanation of this result by expansion of the particle trajectory in 
the random displacement and statistical averaging of a non-perturbative resummation of an infinity of graphs of 
correlations along the trajectory. 

For r] = 0.99 and p ^ 1 one notes a flattening of the scattering function with recovery of the exponent 2/3 power 
law at smaller rigidities. This flattening is definitely present for rj = 1 (no mean component of the magnetic field), and 
correspond to the phenomenological Bohm diffusion regime, as will be seen further below; however, it only extends 
over slightly less than a decade in rigidity for 0.1 ^ p ^ 1, even though for that simulation the dynamic range was 
very large fcmax/fcmin — 10*. At maximum pitch angle scattering, i.e. when p ~ 1 and 77 ~ 1, the scattering function 
g ~ 0.5, i.e. the pitch angle scattering time Tg is of order 2 Larmor times ^l- It should be noted that we define the 
rigidity with respect to the maximum scale of turbulence, which strictly speaking does not coincide with the coherence 
scale /coh of the turbulent magnetic field. In effect, the spatial correlation function of the turbulent component is 
defined as 

(B(x + r,B,x))^(B")i^j^^. (13) 

with S{k) = fc^(B^(fc)) the power spectrum. This integral cannot be integrated analytically for a power-law spectrum 
S{k) (X k^^ but one can check numerically that the maximum of the correlation function occurs at scale Zcoh — 
0.77L,„ax/27r. 
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FIG. 3: Behavior of the averages (Ax^) / At in units of tlc, as a function of the time interval At in units of ti,, for various 
turbulence levels (p — 0.848), and for both the transverse displacement (lower thin line curves) and parallel displacement (upper 
thick curves). One sees the transition from the weakly perturbed propagation regime (Aa;^) oc At^ to the diffusion regime 
(Aa;^) oc At, which appears here as a plateau. The transition duration depends on the turbulence level, and is of order of Ts 
the scattering time. The diffusion coefflcients are given by the levels of the plateaux. Obviously, D\\ S> Di_ for 77 < 1 and the 
two meet in the limit 77 ^ 1, as expected. 



Turnine to the spatial difFusion coefRcients, it is interesting to plot the statistical estimators for D|| and D± given 
by Eqs. (^,(|^) as a function of time for different turbulence levels, and the result is shown in Fig. (||). 

This figure illustrates the transition from the regime in which the particle orbit is weakly perturbed and memory of 
the initial conditions is kept to the regime in which this memory is lost and the particle diffuse, (Ax^) /At ^constant. 
The level of this plateau gives the magnitude of the diffusion coefficient; Fig. also gives an idea of the uncertainty 
in our measurement of diffusion coefficients. Finally, this figure also confirms the expected results ^ D± when 
77 <C I and D\\/D±^ ^ I as 77 ^ I. It should be pointed out that the initial value of the pitch angle cosine was 

/i = I/\/2 in all simulations; we have checked that our results are insensitive to this value as long as the turbulence 
level rj > O.I, as expected. 

In Fig. (^) , we show the behavior of the parallel diffusion coefficient as a function of rigidity for various turbulence 
levels. The dotted lines correspond to the approximation of obtained from the calculation of Tg using Eq. (|), and 
the agreement appears excellent. This study does not confirm the existence of a Bohm scaling. More precisely, the 
Bohm diffusion coefficient Z3b oc ri,v only applies at 77 = I in the range O.I P ^ 1, in agreement with the similar 
conclusion for the scattering function. In all other cases the quasi- linear prediction is verified, i.e. Un cx p^^^ for 
p < 1. We also found that a diffusion regime exists for rigidities greater than the upper bound of the resonance region, 
i.e. p > 1, for as far as we have searched, or about 1.5 decade. In this regime p > 1, Du oc p"^^^, for all values of 77. 



B. The issue of transverse difFusion 



In Fig. (|5|), we plot the behavior of the transverse diffusion coefficient as a function of rigidity for various turbulence 

levels. It is useful to plot also the quantity (^D±/ D^^^^^ as shown in Fig. (^). Indeed, the noise of the simulation is 
then reduced and this figure allows to compare directly the power law behaviors of D± and Z^n . 

This figure indeed reveals a clear trend. For all 77, the ratio D_l/D|| is independent of rigidity for p < I , and 
scales as for p > 1. A similar regime has been found by Giacalone & Jokipii Q for p < I, albeit with slightly 
lower values than ours. This constancy is interpreted in the following as the signature of diffusion due to the chaotic 
wandering of the guide center carrying field lines. The importance of the guiding center diffusion was pointed out by 
Jokipii 1^ as early as 1966 in order to correct the quasi-linear result; however this derivation does not apply to high 
turbulence levels. Finally, the ratio D±/D^\ converges as expected to I for all p when 77 I. However it is interesting 
to note that even at 7; = 0.99, there remains the power law dependence for p > 1, D±/D\\ cx p~^. 

We have found evidence for subdiffusive regimes {Ax^) oc At™, with m < I, at low enough rigidities p < 10^^ 
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FIG. 4: The parallel diffusion coefficient in units of tlc as a function of rigidity for various turbulence levels. The symbols 
and vertical dashed lines are as in Fig. (O). The dotted lines are obtained from the pitch angle scattering rate, using Eq. (0). 
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FIG. 5: The transverse diffusion coefficient as a function of rigidity for various turbulence levels, with the same notations for 
the symbols as in Fig. (^. 

and for rj < 1. On analytical grounds one expects m = 1/2, corresponding to the so-called process of compound 
diffusion and we have found values of m close to this value indeed. However we have not been able to 

investigate in detail this issue, as it is very consuming in terms of computer time. In effect, this can be studied only 
using the GJ algorithm, since it takes place at low rigidities outside the dynamic range of the FFT algorithm. We 
have thus decided to postpone the study of these anomalous regimes to a subsequent publication. 



C. Caracterisation of magnetic chaos 



When the magnetic field is a superposition of a mean field and an irregular component depending on all three spatial 
coordinates, the field line system generically exhibits chaotic solutions. For instance it is sufficient to use a distribution 
of Fourier modes following a power law in wavenumber to obtain a chaotic system. However a two-dimensional field 
cannot have chaotic field lines, and a one-dimensional system cannot produce transverse diffusion, as the particles are 
confined in a flux-tube by conservation of the adiabatic invariant pM . An example of this phenomenon is shown in 
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FIG. 6: The square root of D±/D^^ as a function of rigidity for various values of 77. The notations of symbols are as indicated 
and as in previous figures. The dotted curves overlaid on this figure correspond to the classical scattering result given by 
Eq. (p^), and correspond from bottom to top to the represented values of 77 in increasing order except for 77 = 1. These models 
account marginally for the numerical results for high rigidity particles p ^ 1 and small turbulence levels 77 < 0.5, but diverge 
significantly from the experiment in other regimes. 
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FIG. 7: Transverse displacement of particles in three-dimensional chaotic turbulence (thick line) and in one-dimensional non- 
chaotic turbulence (thin line). In both cases, 77 = 0.5, and for one- dimensional turbulence, the inhomogeneous component 
is taken to depend on the coordinate z, with the homogeneous magnetic field component lying along the z axis. Note the 
difference in behavior; in one-dimensional turbulence, the particle is confined to a fiux tube and does not diffuse. 



Fig. in which we show the transverse wandering of a particle in three-dimensional and one-dimensional turbulence. 

In a three-dimensional chaotic system the separation between two initially adjacent field lines first increases expo- 
nentially oc exp(s//K) as a function of the abscissa s along the field line, with characteristic Lyapunov exponent Ik, 
also called the Kolmogorov length. When the separation has become larger than the coherence length of the magnetic 
field, it behaves diffusively with magnetic diffusion coefficient = ^2As ' ^ti^^e Ar denotes the separation between 
the two field lines. 

Our numerical computation of the field lines clearly displays this two- step behaviour. In Fig. |^, we plotted the 
separation squared between two field lines as a function of the curvilinear abscissa for i] = 0.08. These calculations 
have been obtained by integrating the equations defining the field fines, namely dx/B^ — dy/By — dz/B^, instead 
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FIG. 8: The square of the separation distance between two field fines as a function of the curvifinear abscissa along the field 
line. The exponential divergence followed by the diffusion regime is clearly identified. The transition between these two regimes 
occurs at s --^ imax- 
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FIG. 9: Kolmogorov length and magnetic diffusion coefficient as functions of r). The two lengths are normalized to the largest 
scale of turbulence Lmax- 



of integrating the particle equation of motion. Figure ^ clearly shows this two-step behavior and confirms that the 
transition from one regime to the other occurs when s ^ imax- 

This calculation allows us to measure the two lengths and Dm with a relatively good accuracy. The results are 
reported as function of 77 in Fig. (^). 

The effective transverse diffusion of particles in a chaotic magnetic field has been derived by Rechester & Rosen- 
bluth [p^ . Here we extend their argument by assuming that the primary transverse diffusion is anomalous (sub- or 
super-diffusive). The problem can be stated as follows. After n scattering times, parallel diffusion leads to a diffusion 
in curvilinear abscissa (As^) = 2D\\Tsn, whereas the transverse primary variation causes transverse displacement such 



that (Ax^^) 



^n", with a = 1 for normal diffusion, and a < 1 for subdiffusion. Because of field line exponential 



divergence, until the separation is of order the correlation length, say after scatterings, the transverse displace- 
ment is amplified exponentially by a factor e^''"/'^ with s„ = yJ2D\\Tsn. After Uc scatterings iric ^ 1), an effective 
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transverse diffusion coefficient can then be estimated as: 

(Ax2_) As 



2As At 



(14) 



Because tlie particles almost follow the field lines, the first factor can be approximated by the magnetic diffusion 
coefficient D^, and one gets: 



v/3n72 



(15) 



The number Uc is obtained by equating the separation distance and the transverse correlation length of the field lines 
l± (in our case /_l — ^max): 



tlu^' exp ; 



(16) 



which leads to 



I — 3 , 



(L 



(17) 



where / = vt^ (the scattering length). The main result is that magnetic chaos amplifies the transverse diffusion in 
such a way that it becomes a sizable fraction of the parallel diffusion: 



D, = 



2D„ 







Ik log 







Dn 



(18) 



As can be seen, the primary subdiflfusion does not refrain the effective diffusion due to chaos. When a = 1 (non- 
anomalous primary diffusion), the logarithmic factor reads ~ log {I ±/glK), with g the scattering function as before. 

Finally note that the above regime of diffusion applies at late times after Uc scattering times. The intermediate 
regime, for n scattering times, with n < ric, leads to sub-diffusive motion (compound diffusion), with Ax\ oc At^/^, 
see for instance Ref. [|3|. We have found evidence for such a regime, but a detailed study of its behavior lies beyond 
the present work and is defered to a later study. 

No theory gives the ratio D„^/Ik, except for some toy models such as the Chirikov- Taylor mapping However our 
numerical experiment can provide a fairly accurate estimate of this ratio. In particular we find that the Kolmogorov 



length Ik oc Lmax?? 



-0.9±0.1 



and that the magnetic diffusion coefficient D„ 



T „1.4±0.1 



as long as 77 < 0.5 [see 



Fig.(y)]. Beyond this limit, our calculations of Dm and Ik do not provide accurate estimates of these lengths, 
especially for the Kolmogorov length which loses its physical meaning when 77 reaches unity. Therefore bearing in 
mind that the two diffusion coefficients become the same as 77 — > 1, we conjecture that the result should be: 



n 



2.3±0.2 



Dn 



(19) 



when r/ < 0.5. This non perturbative result would be in agreement with the perturbative result obtained by Chuvilgin 
and Ptuskin Ref. |l^ for small amplitude (written A) large scale varying fields, the ratio between the two coefficients 
being proportional to A*. 

Finally, it is important to note that our numerical results for D±/D\\ shown in Fig. |6|have been obtained indepen- 
dently of the above magnetic diffusion law. We find, in agreement with the above relation, that D±/D^\ is independent 
of p for p < 1, and that D^/D\\ cx rf -^ provides a good fit to the scaling observed at 77 < 1. However the numerical 
prefactor in this relation is rather of order ~ 0.2 for 77 < 1, whereas it should be ~ 1 if the extrapolation could be 
taken up to 77 = 1. Nevertheless, the above provides solid evidence in favor of a dominant contribution of magnetic 
diffusion to the process of transverse diffusion. 



IV. SOME APPLICATIONS 



In this section, we offer revised estimates of the maximal energy that can be attained by Fermi acceleration 
mechanisms by comparing the acceleration time and the time of escape of cosmic rays outside of the accelerating 
region using the results obtained in the previous section. We first consider the case of Galactic supernovae remnants 
(SNR) and so-called superbubbles, and then turn to the case of jets in extragalactic sources. 
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A. Supernovae remnants and super-bubbles 



The lagging questions of the production of cosmic rays in supernovae remnants has been recently reviewed in 
Ref. [ p^ . One of the major problems in accounting for the observational data is that the maximum energy achieved 
by the Fermi process in the SNR shock is well below the so-called " knee" -range 10^^ — 10^^ eV. If one uses the Bohm 
approximation to the diffusion coefficient, there is hope to reach the knee energy with sufficient efficiency to expect 
significant 7-ray emission resulting from tt" decay generated by pp-coUisions. The lack of detection of these gamma- 
rays JlOt or their marginal detection at best, ruined these optimistic assumptions. The Fermi acceleration at a shock 
of velocity Us is characterized by an acceleration timescale tpi — 2D/u1. In most of the shock region, D ~ Z3|| hence 
tpi — Ts/ 01 — ti^/{g0i), where (3s = Us/c. The maximum energy is limited by the age of the supernovae remnants, 
and one thus obtains 

esNR ~ 1.8 X 10"Zg ( ' ( ' ( eV . (20) 



10- V VSOOyry Vl^G 

This result differs from ||l8|| o nly by the factor 3g. This factor is close to unity when the Bohm scaling applies; but, as 



we have found in Section III, in fact g oc p^^^. Strictly speaking this scaling is valid for Kolmogorov turbulence, and 
one expects the turbulent magnetic field downstream to differ from isotropic three-dimensional Kolmogorov, but the 
above scaling serves well for order of magnitude estimates. Moreover p <ti 1 for Larmor radii smaller than turbulence 
correlation length which could be the case even for the most energetic particles. The Bohm approximation thus 
appears very optimistic. 

Super-bubbles correspond to huge cavities created by ~ 100 SNR shock waves built around massive stars associa- 
tions. The size of these regions can be a sizable fraction of the galaxy disk thickness h ~ 120 pc). In effect, a typical 
super-bubble radius can be estimated as |^ 

(io^)''^(i£^)"''C^ (21) 

where L measures the mechanical luminosity of the OB-stars association, uq the particle density of the surrounding 
interstellar medium ~ 1 cm"'^, and ^Myr ~ 30 is the super-bubble lifetime in units of Myr. The bubble plasma is more 
dilute than the interstellar medium by at least two orders of magnitude and thus the Alfven velocity is much greater. 
This density reads ||2^ 

~ 1 V in-2r.m-3r6/35 19/35 -22/35 2/7 

n-SB — 1-0 X 10 cm L^^ Uq t-^^^ , {ZZ) 

where kq is a number of order unity pof . The bubble is traversed by many shock fronts propagating with velocity of 
order or greater than the Alfven velocity; a second order type of Fermi acceleration is thus at work. Its acceleration 
time scale is given hy tF2 ~ (c/Va)^Ts = {c/VAf't'L/ Q- The maximal energy is limited by escape of the particles, which 
is governed by the diffusion across the galaxy disk thickness for the most energetic. A strict lower limit to the time of 
escape t^sc can be obtained by using the parallel diffusion coefficient, since Tosc oc l/D with D the diffusion coefficient. 
Transverse diffusion would improve the confinement time and thus lead to a higher maximal energy; however one 
should then take into account the fact that magnetic lines come out of the galaxy disk and unfold in the halo. Let us 
consider the lower estimate: 

9- (23) 



The maximal energy for acceleration by second order Fermi process in super-bubbles then corresponds to Tosc — tp2i 
and reads 



esB^4x 10i2eV5^(^^) tgf, (24) 

where we used Uq — lcm~^. Therefore the second order Fermi process in super-bubbles might cover the knee range 
with slightly optimistic assumptions, since the magnetic field intensity can easily reach 10/iG in these super-bubbles, 
and ^Myr 30. Moreover, at that maximal energy the rigidity reaches unity and therefore g ~ O.Sry, smaller but close 
to unity . At this point it is useful to recall that the confinement limiting energy of cosmic rays in the galaxy obtained 
by equating the Larmor radius with the thickness h is of order Z x 10^^ eV. 
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B. Extragalactic Jets and Hot Spots 



Extragalactic jets emanating from Active Galactic Nuclei have been considered as possible sources of ultra high 
energy cosmic rays with E > 10^^ eV because the confinement limit in a jet of radius Rj and bulk Lorentz factor T is 

e^,^10-^oVZr(^)(^) (25) 



IG J \lpc 

and for the powerful, strongly collimated and with terminal hot spots FR2 jets, the product BRj is roughly uniform 
in the jet and is estimated as 

where is the mass of the central black hole. Indeed, asymptotically, the magnetic field at the edge of the jet is 
dominated by its toroidal component; therefore the product BRj is governed by the current generated by the central 
engine along the axis. It slightly decreases alongthe jet, because the return current progressively establishes through 
wrapped lines off axis like butterfly wings |2^, p2j. These jets are launched if the magnetic field intensity is close to 
equipartition with the radiation pressure in the central region, i.e., within 10 gravitational radii. This corresponds 
to B lkG(A/*/lO^M0)"i/2 within lOAU (M^/lO^Mg). Thus the performance of the jets as UHE cosmic rays 
accelerators tightly depends on the nature of the central engine. At the base of the jet, B ~ 100 mG for Rj = Ipc 
is a reasonable number. In the hot spots of the FR2 jets such as those of Cygnus A, i? ~ 10^"* G for a region of size 
1 kpc. Thus with T = 10, the confinement condition in jets rules out the possibility of generating ultra-high energy 
cosmic rays of energies larger than ~ lO'^" eV. In the case of FRl jets, which are less powerful, less collimated, and 
without hot spots, the limiting energy is even smaller since the product BRj, although not well known, is very likely 
lower. 

As usual, the escape of the highest energy cosmic rays is governed by diffusion across the jet and 

' - ^ ~ J y (27) 



2L>_L " 2 772-3c2 tL ' 

where we have use our previous result that perpendicular diffusion is governed by magnetic diffusion of the field line, 
and D±^ ^ Q.2rj^-^D\\. One thus finds that indeed most high energy cosmic rays escape before reaching the end of the 
jet, since 

to be compared with a travel time of ~ 1 Myr to travel 300 kpc, the typical length of extragalactic jets. Here as well 
the maximal energy for Fermi acceleration is obtained by equating Tcsc with the acceleration timescale for acceleration 
in shocks moving at speed PsC. This gives 



''^,^<^i- (29) 



With the plausible assumption of Kolmogorov turbulence, g ^ Q.hrjp^/^ , we finally obtain the maximal energy as a 
fraction 0^ of the confinement energy, the estimate being weakly sensitive to the turbulence level: 

w^io-«vft»zr(^)(i)^ (30) 

Centaurus A is a well-known example of active galactic nucleus, actually the closest to us (distance 3.4 Mpc), 
which displays FRl non rclativistic jets moving at speeds ~ 5 x lO'^km/s. The jets have both radio and X-ray band 
synchrotron emission with luminosity Lx — lO'^^erg/s extending over several kpc. They have been studied in detail 
with high resolution interferometry |23j and recently with the Chandra X-ray satellite . The radio knots and the 
X-ray knots are identical in the inner jet. The minimum pressure magnetic field is i?cq — 60/iG and the maximum 
Lorentz factor of the electrons 7max — 8 x 10^. The inner jet has radius Rj ~ 30 pc and a constant opening angle 
of 6°. The product BRj ~ 1.8 x 10~3G.pc is clearly too low to produce ultra-high energy cosmic rays. In any case 
it has been shown that even if Centaurus A could accelerate cosmic rays to the highest energies observed ^ 10^° eV, 
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their transport to Earth, affected by diffusion according to the rules derived in this paper would lead to strong energy 
losses by increased travel distance and anisotropy incompatible with present observations [ p5| . 

In the hot spots, the escape is again dominated by diffusion at high energies, but parallel diffusion is more likely 
unless there is no ordered field. For a hot spot of size i?hs ~ afewkpc and magnetic field intensity B ^ 10^^ G as in 
Cygnus A, the confinement limit is 



(31) 



and the maximum energy achievable with a non-relativistic shock is 

Emax - Psg£c\ (32) 

The most extreme energy that can be obtained is when the turbulence is high enough that no organised field is set 
up in the hot spot, and the shock is midly relativistic Ps ~ 1. But synchrotron emission of hot spots, like those of 
Cygnus A, does not favor this view. Indeed the synchrotron emission by relativistic electrons cuts off in the infrared 
range. Since an electron of Lorentz factor 7 synchrotron radiates around a frequency I'syn — 115(-B_l/10~*G)7^ Hz, 
an observational upper bound on the electron Lorentz factor is 



B 



-1/2 



7max^lOM YF4g) ' ^^^^ 



and the corresponding rigidity is 



Now this maximum electron energy is obtained by the same Fermi acceleration process, limited by synchrotron losses. 
We remind that the characteristic time for a synchrotron radiative process of an electron of energy e = ^rrieC^ is 



where ctt is the Thomson cross-section and /Iq the magnetic permitivitty of vacuum. By equating this loss time with 
the first-order Fermi acceleration timescale, we obtain the maximum Lorentz factor that can be achieved 



B 



10-4G 



-1/2 

(36) 



and the scattering function is calculated for the maximum rigidity pe- Therefore, assuming again Kolmogorov turbu- 
lence with g{pe) — O-^VPe^^ 1 obtain the approximate value of the turbulence level by equating the two expressions 

for 7max: 



0.1 J \10-'^G J Vlkpc/ 

Since (3s > 0.1 is very likely, the required turbulence level is rather low. This, in turn, reduces drastically the maximum 
energy of cosmic ray acceleration, using Eq. (|3^), Eq. (|3l|) and knowing that g ^ rip^~^: 



lO^^eV , (38) 



where we explicited the scaling with (3 the exponent of the power spectrum of magnetic fluctuations; however note 
that the estimate of rj must be changed with /3. For Kolmogorov turbulence, /3 = 5/3 and using the upper bound on 
r] Eq. ( p7| ) above, the prefactor is of order 10~^(/3s/0.1)^'^, and acceleration is not sufficient to account for the highest 
energy cosmic rays by several orders of magnitude. 

This limit cannot be circumvented easily, since it is severely constrained by the cut-off frequency of the electrons 
synchrotron emission. The only parameter that could be modified without affecting this cut-off frequency is the 
turbulence index (3. If one considers Kraichnan turbulence /? = 3/2 instead of Kolmogorov turbulence, the prefactor 



16 



{PsV)^ is changed into {PsV)'^, but rj itself is lowered by a factor 10 due to the modified dependence of g on pe- 
Furthermore if hot spots were to radiate synchrotron emission in X-rays, this would increase j^,^^ by a factor 10 only, 
and would not affect drastically our conclusions. Finally the numbers considered above are consistent with recent 
observations of the Cygnus A hot spots by Chandra , which give an accurate measurement of the magnetic field 
intensity ~ 1.5 x 10~^ G to within a few tens of percents, as obtained by the ratio of the synchrotron-self Compton 
luminosity over the synchrotron luminosity, a value which is furthermore close to the equipartition value if there are 



no protons! These measurements also confirm model dependent estimates proposed in 1986 |27 . Finally, as an aside 



the same reasoning allows to estimate the level of turbulence required to get the X-ray emission in Centaurus A: with 
7™„, ~ 8 X 10^ Pe ~ 1.2 X 10-4, and rj ~ IQ-^Pl 



V. CONCLUSION 



Let us first summarize the results we have obtained. The scattering function g has been found to follow the scaling 
predicted by quasi-linear theory in the inertial range Pmin < p < 1 for weak to strong turbulence. However we found 
that scattering still operates for p < p,nin, contrary to the predicted sudden drop of the scattering function; this 
facilitates the injection of particles in Fermi processes. For Larmor radii larger than the correlation length p > 1, 
scattering decreases as a power- law in rigidity unlike the predicted sudden drop of t;. Therefore high rigidity particles 
still diffuse. One should also mention that the lack of scattering encountered in weak turbulence theory for particles 
having pitch angle close to 90*^ is cured in strong enough turbulence. 

The perpendicular diffusion turns out to be very different from the prediction of the quasi-linear theory. Our 
investigation of the chaos of magnetic field lines characterised by a Kolmogorov length and a diffusion coefficient with 
space increment indicates that this process of magnetic diffusion governs the transverse diffusion of particles. 

Our numerical experiment shows that the phenomenological Bohm approximation, characterized here by g ~ 0.5 
and D = asfLV with ae ~ 0.7, only applies in a limited range of rigidities 0.1 ^ p ^ 1., and only in the case of pure 
turbulence rj — 1. Many estimates in astroparticle physics, that rely on the Bohm conjecture, must be reconsidered. 

The slow decrease g oc p~'^l^ of scattering for cosmic rays with Larmor radius larger than the correlation length of 
the magnetic field, which implies D oc p' 1^ ^ is of potential importance to the transport of high energy cosmic rays in 
our Galaxy as well as ultra-high energy cosmic rays in the intergalactic medium. 

The accurate knowledge of the transport coefficients allows to be more conclusive than before on the performances 
of Fermi acceleration in some astronomical sources of high energy cosmic rays such as supernovae remnants, super- 
bubbles and extragalactic jets. Using new Chandra data, the turbulence level and the maximum energy for electrons 
and for cosmic rays can be determined. We confirm the difficulty to obtain energies larger than 10^^ eV in supernovae 
remnants and shows that the " knee" range of the cosmic ray spectrum could be accounted for by second order Fermi 
acceleration in super bubbles. We also confirm that FRl jets, such as Centaurus A, although radiating synchrotron 
in X-rays, cannot produce UHE Cosmic rays. On the contrary, FR2 jets can produce cosmic rays up to 10^° eV, 
but presumably not more, owing to a fairly good confinement; however most high energy cosmic rays escape before 
reaching the end of the jet. Hot spots of powerful radio-galaxies have always been considered as a promising source, 
but we have found that, because of the low turbulence level implied by the synchrotron cut off frequency, cosmic rays 
escape rapidly along the mean field lines by fast parallel diffusion and acceleration is not effective above emax ~ 10^'* eV 
for a shock velocity /3s ~ 0.1. 

Our paper left opened several important issues that we are currently investigating. In particular it seems crucial 
to investigate in more detail the existence of subdiffusive regimes at low ridigities for which we have found evidence. 
These regimes play a crucial role in the acceleration processes at perpendicular shocks fl^ . Furthermore, we have 
described magnetic turbulence as an ensemble of magneto-static modes distributed according to a power law spectrum. 
This approximation is justified by the small Alfven velocity when compared to the velocity of the particles. However 
it would be interesting to investigate the effect of temporal and spacial intermittency on the transport properties. 
Finally we are currently investigating the transport properties of particles in non-isotropic turbulence as may be 
encountered in the vicinity of a shock wave, in particular in the downstream medium. The consequences on Fermi 
acceleration will be presented in a forthcoming paper. 



APPENDIX A: THEORETICAL APPROACH 



The diffusion resulting from the random variations of the momentum due to the irregular magnetic field can be 
formalized as follows. Energy conservation, and thus p conservation, allows to treat the problem as random rotations 
of the unit vector u such that p — pu: u(t) = io)u(to)- Assuming that the correlation functions of the components 
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of u are integrable over a caracteristic time Tg, the diffusion coefficients are given by 



{ui{t)uj{t + T))dr 



The correlation matrix is derived by making the appropriate average, after solving the stochastic equation: 

u = Q(t)u 



(Al) 



(A2) 



where the gyro-matrix Q{t) = sign(q) J^, ba(t) being the reduced components of the magnetic field expe- 

rienced by the wandering particle and is a 3 x 3 matrix, with components {Ji)jk — ^ijk, where e^jfe is the fully 
anti-symmetric Levi-Civita tensor, and £123 = 1. 

Note that the Ja are generators of a Lie algebra such that 



JaJp = e/3 (g) Gq, - SapId 



(A3) 



where Id represents the identity matrix, = —11^ and JaJp — JpJa = ^ap-yJ^y, where {ea} is the orthonormal basis, 

the orthogonal projector over the plane transverse to the direction a. We have bi — hi, 62 = ^2 and 63 = feo + 
with (b^) = T] and &q = 1 — r/. Moreover h{t) = b[xo -I- p£,{t)], with ^ = u. The time variable is measured in Larmor 
time units, the space variables are reduced to Lmax and wave numbers are accordingly dimensionless and varies from 
1 to l/pm, where = fcmin/fcmax — -^min/-^max- We make the three following assumptions: 

i) The random process becomes stationary beyond the correlation time Tg. 

ii) The random process can be approximated by a Markoff process beyond the integration time Tg. 

iii) The random process h{t) is supposed to be specified. For instance, it is gaussian with a known correlation 
function (b(t).b(i')> = r(t - t'). 

Assumptions i) and ii) allow to calculate the correlation function with an average matrix that describes the relaxation 
of the correlations: 



(A4) 



where Rijlr) — {Rij{t + T,t)). Assumption iii) is not exact, of course; however the numerical experiments provide 
correlation functions that allow to get a good " guess" . Then the theoretical method allows to calculate the solution 
through iterations, starting with a gaussian approximation, and then estimating, if necessary, non-gaussian corrections, 
given a skewness factor. 
The formal solution reads: 



Rit,to) = ( Tcxp 



(A5) 



where the symbol T represents the " time ordering operator" that organises the expansion of the exponential operator 
in products of non commutative operators that are in chronological order. The result can be factorized as the product 
of the unperturbed rotation in the mean field times some relaxation operator: 



R{t,to) = Ro{t-to).(Texp 



dTfi(T) 



(A6) 



with n{t) = RQ^{t - to)5n{t)Ro{t ~ to) (see [|8[ |29[ |o| for technical details) 



1. Quasi-linear approximation rj <^ 1 



For r] small and for a broad magnetic spectrum insuring a short correlation time of the random force compared to 
the scattering time, the quasi linear theory applies This allows to make two approximations. First, the relaxation 
operator can be calculated to the lowest order, the so-called "Bouret approximation" ^ 
a summation of all the " unconnected diagrams" : 



30 , which corresponds to 



jno+M)t 



(A7) 



18 



with 

M= f dTy2(bo.{t)ba{t~T))jMr)Jc. (AS) 
The simplest way to derive this result is to linearize the evolution equation for R(t, to): 

^SR^noSR + SnR+ ... ; (A9) 

which one then solves to lowest order for SR as a function of Sfl and inserts the result in the evolution equation for 
R. Then we use the isotropy of the spectrum to obtain 

M = - dTr(r) V JMt)J^ . (AlO) 

The sum of operators can be simplified to: 

JaRo{T)Ja = — 2coscjoTn|| — (1 + cosa;oT)nj^ — sincjoTJs , (All) 



where ujq is the reduced Larmor pulsation in the mean field, thus = ^/T—r]. 

Second, the correlation function of the magnetic irregularities experienced by the particles is calculated with un- 
perturbed trajectories. 

r(r) ^ ro(r) = I ^53^(k)e^'"'-«o(-) (A12) 

where ^o(''') = /J" dT'i?o (''■') •u(O), and Ssd, as the notation indicates is the three-dimensional power spectrum in 
Fourier space. These two calculations, thanks to commutation properties, lead to a matrix of the form 

R{t) = Ro{t) exp [-.giiHiif - g^n^t - goJst] (A13) 

The factors g are small numbers of order rj that contrains the usual resonances of the quasi linear theory in TT5{knp\i.j,\ — 



nwg). These resonances come from the cosine and sine factor s in E q. (All) and from the expansion into a Fourier 



sequence in moo of the exponential involved in ro(T), see Eq. (A12), which introduces Bessel functions of all orders. 
However, in pratice, only the main resonances for n ~ ±1 are retained because the higher resonances involve shorter 
and shorter wavelengths which contain less and less energy for usual power law spectra. The contribution in J3 
modifies the gyro-pulsation in the rotation matrix Ro{t) and therefore is unimportant. We finally retain the following 
result: 

R{t) = e-9"*n|| + e-<^^'R^{t) , (A14) 

where Rgit) is the product of the rotation and the transverse projector. We thus obtain the unexpected result that 
the transverse relaxation is longer than the parallel relaxation since g|| — 2g±. Then the correlation functions are 
obtained by averaging over u(0) and there comes a major problem of quasi-linear theory because the functions g are 
proportional to ri{p\fi\)^~^ for > fim = Pm/p and vanish for < because of the lack of resonance. This 
introduce long tail contributions to the correlation functions. This is the symptom of the "sticky" regime for pitch 
angles close to 90° that tends to dominate the diffusion coefficients; which requires to take into account mirroring 
effects and/or overlapping of the resonances for /i > close to pm and those for ^ < close to —pm, as suggested 
in pl[|. This difficulty disappears in strong turbulence and for large enough Larmor radii. 



2. Theoretical hints with no mean field 



The fully deductive theory of this regime is quite difficult. However some attempt can be proposed in the case 
where g < 1, i.e. < ^s, when the correlation time [decay time of T{t)] is shorter than the scattering time. 
Thus, for a time longer than the correlation time, we can keep part of the quasi-linear theory, namely the expression 
of the relaxation operator involving the integral over T(t). Technically, this corresponds to the summation of the 
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unconnected diagrams of the expansion of i?(t,io) in Eq- (A5), the other diagrams ("nested" and "crossed") being of 
smaher orders. Therefore the correlation functions (t) asymptotically decay like e~^** and 



9-9* 



/ r(r)dT 
Jo 



(A15) 



Now the main difference comes from the estimation of the correlation function of the field experienced by the particles: 

d^fc 



r(r) 



(27r)^ 



.53D(k)(e''"'-«(^)) 



(A16) 



We propose the following heuristic estimate. Because the particles follow the field lines when their Larmor radius is 
smaller than the wavelength of the modes, we consider only the modes such that kp > 1. The dominant contribution 
in the averaged exponential is then for short time, ^(r) ~ tu(0). Because of the random distribution of u(0) over the 
unitary sphere, we get 



r(r 



) =^ / Sik) 

Jkp>l 



sin kpT 



lkp>l kp'T 

A similar result is obtained with a gaussian evaluation of the average 



dk . 



k'p't / C{t)At 







Inserted into the integral over the spectrum (restricted to fc > l/p), it leads to: 



r(r) 



k>l/p 



dkS{k) exp 



ifc^V fc{r')AT' 



(A17) 



(A18) 



(A19) 



In the integral, C(r')dr' can be approximated by t for r < Tg (= \/ g in reduced units) and by 1/g for r > Tg. 
Therefore 



AkS{k) 



k>l/p 



VStt ^ f kp 
2kp 



3g 



e 



(A20) 



where $(x) — e ^ dy. When g is small, because kp > 1, we get a simple result close to the previous one: 



d.^ 

k>l/p '^P 



(A21) 



This gaussian evaluation indicates the error made by the previous assumption. Thus, for small p, we obtain the 
extension of the quasi- linear result, namely g ~ and for p > 1 to g ^ 1/ p. These two approximations are in 

agreement with the numerical experiments, except that the measured drop is in p^^'^ instead of p~^ here. Now the 
range of p values where g is on the order of unity corresponds to the "Bohm estimate" , which is, in fact, the maximum 
value of g achieved for p ^ 1 only. We thus propose the following final estimate for the scattering function: 



kp> 



S{k) 

1 kp 



dk . 



(A22) 



The error on the coefficient is of order ten percent. 
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